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Abstract. In recent decades, there have been many attempts to construct symplectic 
integrators with variable time steps, with rather disappointing results. In this paper we identify 
the causes for this lack of performance, and find that they fall into two categories. In the 
first, the time step is considered a function of time alone, A = A(f ) . In this case, backwards 
error analysis shows that while the algorithms remain symplectic, parametric instabilities arise 
because of resonance between oscillations of A(r) and the orbital motion. In the second category 
the time step is a function of phase space variables A = A(q,p). In this case, the system of 
equations to be solved is analyzed by introducing a new time variable X with dt = A(q,p)dt. 
The transformed equations are no longer in Hamiltonian form, and thus are not guaranteed 
to be stable even when integrated using a method which is symplectic for constant A. We 
analyze two methods for integrating the transformed equations which do, however, preserve the 
structure of the original equations. The first is an extended phase space method, which has been 
successfully used in previous studies of adaptive time step symplectic integrators. The second, 
novel, method is based on a non-canonical mixed-variable generating function. Numerical trials 
for both of these methods show good results, without parametric instabilities or spurious growth 
or damping. It is then shown how to adapt the time step to an error estimate found by backward 
error analysis, in order to optimize the time-stepping scheme. Numerical results are obtained 
using this formulation and compared with other time-stepping schemes for the extended phase 
space symplectic method. 
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1. Introduction 

Adaptive time integration schemes for ODEs are well established and perform extremely well 
for many applications. However, for applications involving the integration of Hamiltonian 
systems, there are good reasons for using symplectic integrators[l]. This is particularly true in 
applications such as accelerators, in which very long orbits must be integrated. Symplectic 
integrators are also used in particle-in-cell (PIC) codes for plasma applications, where for 
economy integrators with lower order of accuracy must be used, but it is important that the 
Hamiltonian phase space structure be preserved. The same is true for molecular dynamics 
codes[2, 3], which are of use for dense plasmas as well as other applications. Similarly, for 
tracing magnetic field lines (e.g., in tokamaks or reversed field pinches) or doing RF (radio 
frequency) ray tracing, symplectic integrators can be very useful. To our knowledge, the 
symplectic integrators in use for plasma and particle beam applications all have uniform time- 
stepping. However, for these applications there are situations in which variable or adaptive 
time stepping could increase efficiency greatly. For example, in accelerators, particles can 
experience fields that vary over a wide range of length scales. The same is true for for PIC and 
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MD codes, for example for particles moving in and out of shocks or magnetic reconnection 
layers, or undergoing close collisions. 

There have been studies of symplectic integrators with variable time steps, but the early 
results were not promising [4, 5, 6, 7, 8, 9, 10]. Among these studies, there have been two 
types of variation of time steps. In the first set of studies, the time step varies with time 
explicitly, and in some cases varies from At = Aj to Af = A2 on alternate blocks of time steps. 
In these studies problems were observed to arise. In the second set of studies, the time step was 
chosen with Af = A(q,p), where q,p are the dynamical variables. For any adaptive scheme 
for an autonomous system of equations, the time step is of this form. For the A(q,p) case, 
the equations are no longer in canonical Hamiltonian form (but may be Hamiltonian in non- 
canonical variables; see below) and indeed unfavorable results are found. The disappointing 
results for both of these cases have contributed to the general impression that if you need 
an adaptive time step integrator, you are better served using a high order non-symplectic 
integration method, f 

In this paper we take a new look at symplectic integrators with adaptive time stepping. 
Our focus is not on developing the most efficient or accurate symplectic integrators, but to 
understand and solve the problems that have been encountered in using symplectic integrators 
with variable time steps, and to outline a method of obtaining an optimal time-step adaptation 
scheme. 

In Sec. 2 we start by considering time steps depending explicitly on time alone as 
At = A(t) = Ao(l +ecosa)f). In Sec. 2.1 we consider several examples of first and second 
order symplectic integrators applied to the harmonic oscillator, and apply modified equation 
analysis or backward error analysis[\A, 15, 11, 12]. This allows us to find the modified 
Hamiltonian system which the integration scheme, including the effect of numerical errors, 
actually solves. From this analysis we find that resonances between co and the oscillator 
frequency Q.q drive parametric instabilities. The resonances have CO = mClo, the lowest order 
resonance having m — 2. For first order schemes the resonance width scales as eAo, while 
second order schemes have width °< sAq. These parametric instabilities explain the problematic 
results obtained in the papers dealing with A = A(f)[5, 7, 8, 9]. In the context of the harmonic 
oscillator, these problems were identified as arising from parametric instabilities by Piche[16]. 

We continue in Sec. 2.2 by numerically investigating time step variations A(f) = 
Ao(l +ecoso)f), by integrating the harmonic oscillator with a symplectic method. We 
find parametric instabilities with resonance widths as predicted by the modified equation 
analysis. We show results of numerical integrations of a cubic oscillator with the nonlinear 
potential V (q) = q 2 /2 + q 3 /3, in which the frequency varies with respect to the action variable, 
Q. = £l(J). The parametric instabilities seen globally for the harmonic oscillator show up 
as nonlinear resonances or islands, where co — mQ.(J), with m depending on the integration 
method and the Hamiltonian being integrated. We observed m = 1 islands for the cubic 
oscillator, integrated with the Crank-Nicolson (CN) method. The calculations of this section 
illustrate the point that, when considering A = A(f ), problems arise because the equations being 
solved become unstable (but still Hamiltonian), and not because the integration scheme fails to 
be symplectic. 

In Sec. 3 we study time step variations A = A(q,p). With A(q,p) °< l/p(q,p), the 
straightforward substitution p(q,p)dt = d% leads to equations (having T as the independent 

f An exception to this is the extended phase space approach of Hairer[l 1] and Reich[12], which we discuss in later 
sections. Some authors[13] also suggest that methods which are reversible under some symmetry R (RT(h)R = T(h)^ 1 , 
where the one-time-step map is T(h) and R is an involution, R = the identity) should be considered in place of 
symplectic integrators. We do not consider reversible methods here because their advantages are limited to systems 
where the original flow also has this symmetry property. 
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variable) that are not Hamiltonian equations in canonical form. Symplectic integrators therefore 
should not be expected to give stable results in this case, and again it has been observed that 
they do not. In Sec. 3.1 we review one approach to solving this problem, by extending the 
phase space (q,p) —> (q,qo,p,po), where qo — t and po is its conjugate momentum. After this 
transformation to an extended phase space, it is easily shown that the resulting equations are 
Hamiltonian in canonical form and can be integrated by any fixed time step (At = h = const.) 
symplectic integrator in this extended phase space. This is a well-known method, and has 
been used to obtain symplectic integration schemes with variable time steps by Hairer[l 1], and 
also Reich[12]. Hairer described this approach as a "meta-algorithm" because any symplectic 
integrator can be used in the extended phase space. We then present numerical results in Sec. 3.2 
for the cubic oscillator (with two different symplectic integrators) which show that, even with 
a step size variation that appears to have an m — 1 resonance, evidence of non-symplectic 
behavior occurs and no resonant islands form. 

In Sec. 3.3 we introduce a second method for dealing with the problems that arise when 
A = A(q,p). The point of introducing a second method is to illustrate that the main problem 
with A(q,p) is that the equations are no longer Hamiltonian, and that there are potentially 
many ways to address this problem. The extended phase space method "fixes" the equations 
by embedding them in a higher dimensional Hamiltonian system. This alternative approach 
does not extend the dimension of the phase space, but rather recognizes the fact that, with z 
as the independent variable and p(q,p)dt — dx, the equations are Hamiltonian equations in 
non-canonical variables[\l\. (This is true only for one degree of freedom, i. e. 2D phase 
space.) That is, the equations can be expressed in terms of a non-canonical bracket[17]. We 
write the non-canonical variables as (x,y) with the numerical scheme giving the mapping 
(x,y) = (x(t),v(t)) — > (X,Y) = (x(x + h) ,y(x + h)) . Equivalent to the non-canonical bracket, 
the fundamental two-form is p (x,y)dx Ady rather than the canonical form dqAdp. We describe 
a method for constructing a non-canonical generating function F(x,Y) which gives a map 
7> (/i) which is a first order accurate integrator that preserves the two-form, and discuss its 
relation with transforming to canonical variables (guaranteed to be possible by the Darboux 
theorem[17].) We then find a complementary non-canonical generating function G(X,y) which 
leads to another non-canonical map Tq{H), which is also first order accurate and preserves 
the two-form. Maps preserving the two-form p(x,y)dx Ady are called Poisson integrators 
[18, 19, 20, 21] and conserve the measure / p dxdy rather than the area / dqdp. We show that 
T F {h)- 1 = T G {-h) and vice-versa, so that the composition T(h) = Tc(h/2) o Tp(h/2) is time 
symmetric (T(h)~ l = T(h)) and therefore second order accurate. Finally, we show numerical 
results in Sec. 3.4 using the composed map T(h), applied to the cubic oscillator. As with 
the extended phase space method, this method is stable, with errors in the energy remaining 
bounded. 

In Sec. 4 we study the choice of time step variation p(q,p). We begin by using backward 
error analysis for a slightly different purpose than in Sec. 2, namely to estimate the form of the 
error of an integrator, and then devise a scheme to minimize the total error. This scheme takes 
the form of an equidistribution principle, of the form p opt (q, p)dt = dx, with uniform steps 
in x. We then compare the error in numerical integration using p op , to the error obtained for 
fixed time steps, to the error obtained for time steps designed to give equal arc length per step, 
and to time-stepping using a weighted average of p upt and p = 1 . We show numerical results 
showing that, among these options for time step variation, p opt does indeed minimize the error. 

In Sec. 5 we summarize our work. In Appendix A we consider the relationship between 
the cases with A(t ) and A(q,p). We show that, while a simple argument motivates the use of 
A(f ), this is only valid for short times. We explain how this is analogous to the problems which 
arise in perturbation theory, where a perturbation causes a frequency shift, but at lowest order 
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it appears a secularity or an instability. In Appendix B we review non-canonical variables, 
brackets, and forms, material necessary for Sec. 3.3. 



2. Integration with variable time step A(r) 

In this section we will investigate symplectic integrators with a variable time step A = A(f ), in 
order to understand better the results that have been obtained in the literature. Using modified 
equation analysis, we will see that the problems that have been reported are best understood as a 
parametric instability due to resonance between the system dynamics and the time dependence 
ofA(f). 



2.1. Modified equation analysis and parametric instabilities 

We first look at the Harmonic oscillator H = p 2 /2 + Q^q 2 /2. With the change of variables 
p — > pj x/Ho an d q —> V&oq we obtain H — £lo(p 2 + q 2 ) /2. The leapfrog (LF) method, which 
is first order accurate when q and p are not taken to be staggered), with A = A(f ) gives 

=qi + A(ti)D.opi, (1) 

Pi+\ = Pi - A(t,-)£2o?i+i- (2) 

Modified equation analysis begins by expanding q, + i = q\ + A(f,)<j , ! - + A(f,) 2 ^,/2 and 

similarly for pt+i. We obtain 

A(t)£ll , i, 

q = ^j-^q + £l p + O(A 2 ), (3) 

p = -Claq - -^Y^P + 0(A 2 ). (4) 

(Terms proportional to A appear at the next order for first order leapfrog. For other integration 
schemes, the A terms will similarly appear at one order higher than the order of the method.) 
This system of equations arises from the Hamiltonian 

„> Q / 2 , 2\ , A(Q^Q 

As expected, the modified equations are in Hamiltonian form [22, 23, 15, 24, 25, 12]. 
With the action-angle transformation q = \/2Jcos(j), p = — \/2Jsm(j> and setting A(f) = 
Ao(l+£cos©f) we obtain 

H' = Q. J ^1 - sin 20^ - n °^ Q£ 7 [sin(20 - cot) + sin(2</> + 0)t)] .(6) 

Keeping only the resonant O(Ao) term °< sin(20 — cot), we obtain 

H' « Cl J - Op-/ °^° £ sin(20 - cot) . (7) 

The canonical transformation y/ = — cot/2, K' = H' — coJ/2 (J unchanged) transforms to a 
rotating frame in phase space and gives 

, / co\ £2nAn£ 
K'= (D.0- -Jy--5_ u _y sin 2 V /. (8) 

We conclude from this backward error analysis that the harmonic oscillator with time step 
A(f) = Ao(l +£ cos oat) is unstable if 

£^A £ 



2 



< (9) 
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since K' is hyperbolic in this case. This instability is a parametric instability associated with 
the resonance CO — 2Qq. Figure 1 shows this predicted region of instability as a function of 
CO and e. Numerical results for comparison are shown in the next section. Carrying out this 
analysis to one higher order in Ao gives an 0{A 2 ) ) correction to the 0- and time-independent 
terms in K', which gives rise to a positive 0{A 2 ) ) shift in the resonant frequency independent of 
e, discussed in the next section. 

Another integration scheme is symmetrized leapfrog, which is second order accurate, and 
takes the form 

A(t/)flo 
q* = «H 2 — Pu 

Pi+i=Pi-A(ti)Q. q*, (10) 
A[f,-)i2 

This leads to a modified equation 



qi+i = q*-\ r — Pi+U 



A(t) 2 Q.l A(f) 2 £2^ 
q = a oP -^i-^p, p = -n q~^r^q, (ID 



with modified Hamiltonian 



12 

A similar analysis to that above for first order LF leads to 

, / Q? n Al co\ fl^£ 
K'=Ul + ^J>--)j+^^Jc OS 2 W . (13) 

Note that for this method, the frequency shift enters at the same order as the resonant term. We 
conclude that parametric instability occurs for 



i^Ag CO 



CllAle 

< ° o ° . (14) 



24 2 } 

This instability is centered around CO — m(i2o +Q.qAq/24), again with m = 2, but in a narrower 
range ~ gA^. (Resonances with other values of m enter at higher order in Ao and e.) Also 
note that, because the frequency shift enters at the same order as the resonant term, the shift is 
comparable to the resonance width for e ~ 1 . 

As another second order accurate example, we consider Crank-Nicolson (CN), 

qi+i = qi + n{ti)Llo , p i+ i=pi-h(ti)klo — ^ — ' ^ 5 ^ 

for which the modified equation analysis gives 

q = ( 1 - ^tV^J P = ~ i 1 - J &oq (16) 

or 

"'K 1 - 4 ^)^^)'* ,17) 

Here, H' is independent of and the action J — (p 2 + q 2 )/2 is an invariant to order Aq, and 
the errors are all in the phase (j>: the frequency Q.(t ) = 1 — A(t ) 2 Q.q/12 is downshifted for this 
case. We conclude that parametric instabilities cannot occur to this order for CN. It is easily 
seen that CN preserves the action J exactly for the harmonic oscillator; this shows that we have 
H 1 = Q.(t)(p 2 +q 2 )/2 to all orders. Further, the explicit time dependence disappears if we 
introduce T such that Q.(t)dt = dr. 
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2.2. Numerical results with A = A(f ) 

The parametric resonances described for the harmonic oscillator in the previous section were 
studied numerically, and the results are shown in Figs. 1 and 2. Figure 1 shows the numerical 
results, using the first order LF method, again with time steps given by A(f) = Ao(l +£cos cot). 
The red lines show the boundaries of the unstable region in co, £ parameter space, as predicted 
in Eq. (9). The colors show the logarithm of the value of the Hamiltonian at the end (f = 300) 
of a numerically computed orbit. For unstable parameter values, the orbit spirals outward, and 
the final value of the Hamiltonian becomes exponentially large. The numerically computed 
region of instability agrees well with the prediction of Eq. (9), and also shows the 0(Aq) offset 
in frequency. In Fig. 2 we examine how, for parameters where the integration is unstable, a 




e 



Figure 1. Numerical results showing a region of instability for a range of values of the 
parameters e and CO, for first order LF. The red lines show the predicted stability limits given in 
Eq. (9). The colors show the logarithm of the Hamiltonian function H evaluated on the orbit 
after integrating the equations to t = 300. This value should be approximately constant along 
the orbits, but it becomes large on orbits that spiral out to large values of q and p. Parameters 
used: Q.q = 5, Aq = 0.05. Initial conditions: q = 0.5, p — 0. The small upward shift in CO comes 
in at second order in Aq. 



circle of initial conditions is stretched out into a long thin ellipse. Since the LF method is still 
symplectic, the area of the ellipses is preserved. 

As described above, first order LF integration of the harmonic oscillator with variable 
time step A(f) = Ao(l + ecos cot) was found to be unstable near the resonance CO — 2Q.Q. For a 
nonlinear oscillator, the oscillation frequency Q.q depends on the amplitude. We thus expect 
that the numerical integration of a nonlinear oscillator would exhibit resonant islands, rather 
than exhibiting global parametric instability. 

For example, consider a Hamiltonian with a cubic potential 

H(a,p) = l±f + t. (18) 

This Hamiltonian has an O-point at (q,p) — (0,0) and an X-point at (—1,0). Inside the 
separatrix of the X-point, this system exhibits nonlinear oscillations. For this system and 
first order LF, the time step variation A(f) = Aq(1 + ecos cot) introduces an m = 1 resonance 
(co = Q.q), which is readily understood by backward error analysis. Specifically, this analysis 
gives a term in the modified Hamiltonian which is proportional to sin <p cos cot, which leads to 
the m = 1 resonance. 
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Figure 2. Numerical results illustrating that, while the first order LF method is area preserving, 
variable time steps can cause a parametric instability. The black ellipses are generated by taking 
a set of points, and propagating them forward in time using first order LF. The black dots show 
one of these points, which is plotted for various times (as labeled). The grey curve shows the 
coordinates (q(t),p(t)) of this point, for each of the computed time steps. Parameters: e = 0.4, 
m = 10, 12 = 5, A = 0.05. 




-1 -0.5 0.5 

q 



Figure 3. Surface of section plot for the cubic Hamiltonian of Eq. (18), with variable time step 
A(f) = Ao(l + ecos £0/). The resonance with the time-step variation causes an m = 1 island 
to open. Parameters used: Ao = 0.3, £ = 0.6, CO = 0.95. The dots are points in the surface of 
section (computed using CN), and the thin line is the separatrix of the X-point at ( — 1 , 0) . 



Figure 3 shows the results of integrating this system with the CN method, and the same 
time steps variation. The frequency CO, < CO < 1, was chosen so that it would resonate with 
an orbit inside the separatrix. (The frequency varies with amplitude from the O-point from 
a maximum £2 = 1 at the O-point to £2 = at the separatrix.) The dots in Fig. 3 show the 
Poincare surface of section for a time-T map, with T = 2n/co. The elliptic point (O-point) of 
a resonant m — 1 island is clearly seen near (q,p) = (—0.35,0). We see that indeed the time 
dependent step size variations lead to an unphysical parametric resonance for the Hamiltonian 
system in Eq. (18). 
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3. Integration with variable time step A(q,p) 

We now turn from considering the time step variation A(f) to variations of the form A(q,p). 
This form of variation is perhaps more natural, since it arises when trying to optimize the time 
step to reduce errors in the integration method applied to an autonomous system. (We will 
return to this issue in Sec. 4.) In this section, we first discuss the problems that arise with this 
variation, and then we describe two methods by which these problems can be resolved. 

Previous attempts to apply this type of time step variation to symplectic integrators have 
suffered from one main, although sometimes unrecognized, problem: the equations are no 
longer Hamiltonian. As has been pointed out[l 1, 12], integration of Hamiltonian's equations 

with variable time steps satisfying p(q,p)dt — dx is equivalent to changing the time variable t 
to x, and integrating the following equations with fixed size steps in t; At = Ir. 

dz^ _ J_ dH(z,t) 
dx p(z) dzj 
dt 1 



dx p(z)' 

where £12 = —£21 = 1, £11 = £22 = (summation assumed). Here, z = (q,p). As they stand, 
these equations are not in canonical Hamiltonian form. Because of this, the application of a 
symplectic integrator should not be expected to give results that are symplectic, since Eq. (20) 
is not a canonical Hamiltonian flow. 

Once this problem has been identified, there are various approaches that can be used to 
find an integrator which preserves the Hamiltonian nature of Eq. (19), even with time steps 
given by p (q, p)dt — dx. In the remainder of this section, we review and discuss two methods 
for integrating Eqs. (20) in such a way that the symplectic nature of the original equations 
[Eq. (19)] is preserved. The first method is to embed these equations into a larger phase space, 
considering t to be another coordinate. This extended phase space approach goes back to 
the work of Sundman [26], and has more recently been applied to numerical methods, both 
non-symplectic [27] and symplectic [11, 12]. In the extended phase space, the new equations 
can again be written in canonical Hamiltonian form, and thus integrated with any (fixed x 
step) symplectic integrator. This extended phase space method is not the only way to integrate 
Eqs. (20) symplectically, however. A second method which we consider in this section (which 
only applies to one degree of freedom problems) is to recognize the the equations can be 
written in terms of a non-canonical bracket. Any method which preserves this bracket will then 
have all the nice properties of a symplectic integrator for a canonical system. We present a 
novel generating function approach for constructing such a method. 

For the numerical results presented in this section, we will consider the Hamiltonian with 
a cubic potential given in Eq. (18), integrated with the two different methods. The variable 
time steps A = A(q, p) are chosen so they might appear to have an m — 1 resonance, since that 
is what is expected to give islands similar to the previous results of Fig. 3. Let 



A(q,p) = A()(l+aiq + a 2 p) . (21) 
For the simulations reported in this paper, a\ = 0.5 and a 2 = 0.25.f For p(q,p)A(q,p) = h, 

t We take ai ^ in order to make the equations asymmetric in p. Otherwise, the symmetry of CN makes the results 
appear symplectic, when in general they would not be. 
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i.e. At = A(q,p) and At = h, we have 

P{QtP) — ~tt~ c = t n , h , y ( 22 ) 

A(q,p) A (1 +a\q + a2p) 

The step size At = h was held fixed, and Ao was adjusted so that after a given number N 

periods of the orbit, t(N) equals x(N). This is done so that the differences between the original 

and extended phase space results are due to the step size variation, rather than a change in 

the average step size. Integrating the original Hamilton equations (19) with fixed step size 

CN (i.e., with a\ = a^ = 0) gives a result with bounded errors in energy [see Fig. 4], as is 

expected for symplectic integrators. However, integrating Eqs. (19) using CN with step size 

given by Eq. (3.2) gives errors in energy with linear growth, as also seen in Fig. 4. This lack 

of boundedness of the energy is to be expected, since varying the step size in this way is 

equivalent to integrating Eqs. (20) with fixed step size, and (20) is not a Hamiltonian system. 

Thus there is no reason that using a symplectic method such as CN would give results with 

bounded errors in energy. (Note that the Hamiltonian is just a useful diagnostic; the underlying 

problem is that phase space area is not conserved.) 



x IO- 



cd 
I 




Figure 4. CN integration of the cubic Hamiltonian [Eqs. (18) and (19)] is shown in gray, and 
the time transformed version of the cubic Hamiltonian [Eqs. (18) and (20), with p as given in 
Eq. (3.2)] is shown in black. Parameters: h = 0.05, A = 1.05134, a { = 0.5, a 2 = 0.25. Initial 
conditions: (qo,Po) — (0.3,0). 



3.1. Extended phase space method 

By a well-known procedure [11, 12, 26, 27], we introduce qo~t and its canonically conjugate 
momentum po and consider the Hamiltonian in extended phase space: 

K(q,p,q ,po) = -j^{H(q,p,q )+p ), (23) 

again with z = (q,p). The canonical equations of motion are 

dq 1 dH(q,p,q ) , , d ( 1 

Tt =W) dp + {H + PQ) dp\-p 

dp 1 dH(q, P ,q )_ {H+po) d A (24) 



dx p(z) dq dq \p 

dqo dK 1 

dx dp p(z)' 

d Po _ dK 1 dH 

dx dqo p{z)dq Q ' 
The last of these equations implies that if we choose po(x = 0) = — H(q(Q),p(Q),Q), then 
dpo/dx = —dH/dx and H + po remains zero. In this case, Eqs. (24) are identical to Eqs. (20). 



Symplectic integrators with adaptive time steps 



10 



So, with this choice of initial condition for po, the actual orbit is followed and K = 0. Note 
that, for this value of po, the other contours of K = const, are not orbits of the original system. 

With the equations of motion in the form (24), we can then use a canonical symplectic 
integrator such as Modified LF (ML) (a semi-implicit, symmetrized version of LF; see [28, 29]) 
or CN with fixed time step At = h. This is the method used in Refs. [11, 12] to integrate 
Hamiltonian equations of motion with variable time steps for the special case dH/dt = 0, for 
which H and po are separately constants of motion. It should be noted that in a symplectic 
integration of Eqs. (24), H + po will be only approximately invariant, so that the results of 
numerically integrating Eqs. (24) will have small errors relative to the numerical results of 
integrating Eqs. (20). 

Note that if H is independent of t , then one needs only to integrate the equations for q and 
p, since po is exactly constant, and qo = t can be found by post-processing. So, even though 
this method "extends" phase space, the dimension of the system effectively remains unchanged. 
Therefore, there is no possibility of parametric instability if the time step is small enough. 

3.2. Numerical results: extended phase space method 

In this section we report numerical results using the extended phase space method [Eq. (24)] 
to integrate the equations from the cubic Hamiltonian [Eq. (18)], using p(q,p) as given in 
Eq. (3.2). In Fig. 5 the energy errors are compared for (a) CN integration and (b) ML integration. 
In both cases, the black curve gives the extended phase space results, while the gray curve 
gives the result obtained by applying the same method with fixed step size to the original 
Hamiltonian equations. As for Fig. 4, the T step h was chosen to make i(N) equal to t (N) 
after N periods of the orbit. We conclude that there is no secular change in the value of the 
Hamiltonian H(q(t),p(t)). The time step in Eq. was chosen because it appears at first to have 
to potential for an m = 1 resonance; in Sec. 4 we will describe a method to minimize the 
integration errors. 

In Fig. 6 is shown a surface of section plot for the extended phase space CN integration. 
In contrast to Fig. 3, there is no resonant island caused by the time step variation. This is 
consistent with the comment at the end of Sec. 3.1. 

As a final check of the extended phase space method, a comparison was made between 
integrating Eqs. (20) where the step size is determined by the position in phase space as 
A(q,p) « l/p(q,p), and integrating Eq. (19) where the step size is a function of time as 
determined by a reference orbit A(f) °< 1 / p (q re f(t) , p re f(t)) . The reference orbit was obtained 
by solving the extended phase space equations [Eqs. (20)], for some initial conditions (qo,po) 
for which the period is T(qo,po). Time-T surface of section plots for several cases give a results 
qualitatively similar to Fig. 3, with the time dependent step size A(f) = A(q re f(t),p re f(t)) 
giving a resonant island centered at the location of the reference orbit. This behavior is similar 
to the secular and unstable cases of perturbation theories discussed in Appendix A. 

3.3. Non-canonical symmetrized leapfrog 

In this section we discuss a second method which can be used to integrate Eqs. (20), which 
works for one degree of freedom systems (2D phase space). In this section, we change the 
notation for a point in phase space from (q,p) toz— (x,y). This is done to emphasize the fact 
that Eqs. (20) are not in canonical Hamiltonian form. 

Equations (20) are of the form u = m/p, where u is the velocity, with V • m — V ■ pu = 0. 
The measure J v ^ pdV is preserved over any phase space volume V (t ) advected with the flow. 
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Figure 5. Extended Phase Space Method. Comparison of (a) CN and (b) ML. Shown in black is 
the error in energy H(q(t),p(t)) for the two methods. For comparison, in gray is shown the error 
in energy for the same method, but applied to the original uniform time step discretization of 
the Hamiltonian equations. Parameters: h = 0.05, Ao = 1.05117, (qo,Po) = (0.3,0), a\ = 0.5, 
a 2 = 0.25. 
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Figure 6. Surface of section plot, strobed at an arbitrary period T, for the extended phase space 
method (computed using CN). Parameters: a\ = 0.5, a 2 = 0.25, Ao = 0.3. Expansion of the 
scale for results from longer runs show no islands anywhere inside the separatrix for any period 
T. 
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Equivalently, for the time-f map z(0) — > z(f ) we havef 

This measure preservation suggests that a measure preserving integrator applied to this 
system would have the same nice properties of a symplectic integrator applied to a canonical 
Hamiltonian system. 

The equations of motion in the form in Eq. (20) are not canonical Hamiltonian equations. 
However, in one degree of freedom, Eq. (20) is a Hamiltonian system in non-canonical variables 
and can be written in the form 

d f = [ Zi ,H] z , where [f,g] z = * J£|* (26) 
at p(z) ozidzj 

is a non-canonical bracket. That is, it is antisymmetric and direct calculation shows that it 
satisfies the Jacobi identity [/, [g,A]J z + [g, [h,f] z ] z + [h, [f,g] z ] z = 0. 

A method of constructing a symplectic non-canonical integrator in terms of generating 
functions [30] starts with the condition that must hold for the time-f map, which is that it must 
preserve the two-form 

p(x,y)dx/\dy = p(X,Y)dX AdY, (27) 

See Appendix B for more details regarding this condition. First, we define 

rj(x,y)= P p(x,y')dy' (28) 



and 

Z(x,y)= J* pix'^dx'. (29) 

With these definitions, we can write 

p(x,y)dxAdy =dxAdi){x,y), (30) 

p{X,Y)dXAdY =d%{X,Y) AdY (31) 

and therefore Eq. (27) becomes 

dxAdr](x,y) =d$(X,Y)AdY, (32) 

We can consider (X,Y) to be functions of (x,y), and then define p = r)(x,y) and Q = %(X,Y). 
These definitions, together with Eq. (32), imply that the one-form 

co = pdx + QdY (33) 

is exact, i.e. dco = 0. This in turn implies 

, x dF(x,Y) dFtx.Y) 
co = dF(x, Y) = — \^-dx + — k^-dY, (34) 
ox oY 

at least locally. We are led to the relations 

n (x,y) = P = dF( f Y l, (35) 
ox 

^X,Y)=Q= d ^l. (36) 

Given the generating function F(x,Y), we solve Eq. (35) for Y(x,y) and then substitute into 
Eq. (36) to obtain X(x,y). 

t For tracing magnetic field lines by dxj/dC, = Bj/Br (i = 1, 2) the same considerations lead to conservation of flux 
J A B^dA at two ends of a flux tube, and det {dx t {Q /dxj(0)) = B ? (x(0))/B ? (x(Q). 
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In essence, what is going on is this: we are transforming to canonical variables (x,p) and 
(Q,Y), where p = rj (x,y) and Q — S, (X, Y). Then we are performing a canonical transformation 
(x,p) —> (Q,Y) in terms of a generating function F(x,Y), with Q = dF/dY and p = dF/dx, 
and expressing the results in terms of (x.y), and (X,Y). 

The identity transformation has X = x and Y — y so that Eqs. (35), (36) take the form 

dF (x,y) 

n{x,y) = — 3- — , (37) 
dF Q (x,y) 

${x,y) = — j — • (38) 

These relations ensure that drj(x,y)/dy — (x,y) jdx — d 2 Fo/dxdy; by Eqs. (28), (29), these 
both equal p(x,y), which is nonzero. This nondegeneracy is necessary for the inversions in 
Eqs. (35), (36) to work. 

We now wish to construct an 0(h) approximation to F(x,Y) in order to find a time-/i map 
(h = At) for the system given by Eq. (20). Let us try the generating function 

F(x,Y)=F (x,Y)+hH(x,Y). (39) 

From Eqs. (35) and (36) we find 

ox 

t;(X,Y) = ^^±+hu(x,Y), (41) 

where u(x,y] — dH(x,y)/dy and v(x,y) = —dH(x,y) jdx. Substituting Eqs. (37) and (38) into 
these equations gives 



r](x,y) = r,(x,Y)-hv(x,Y), %(X,Y) =%(x,Y) +hu(x,Y). 



(42) 



This is the form of our integrator: we integrate p to give r) (x,y) and | (x,y) as in Eqs. (28) and 
(29), iterate the first of (42) to solve for Y and, after substituting for Y, iterate the second of 
(42) to solve for X. We do not need the explicit form for Fo or the Hamiltonian H. 

In order to see that Eqs. (42) are a first order accurate integrator for the system in Eq. (20), 
we expand to first order in h, and find 

0=(Y-y)^p^-hv(x,y) + O(h 2 ) (43) 
dy 



d^y) 

dx 

With Eqs. (28), (29), this leads to 



{X-x) '" 9 Y' ,jrj =hu(x,y) + 0(h 2 ). (44) 



X=x+h U( ^\+0(h 2 ), Y^y + h^ + Oih 2 ). (45) 

p(*,y) p(*,y) 

This implies that Eqs. (42) do constitute a first order accurate integrator for Eqs. (20). It is 
clear that the scheme in Eqs. (42) is the non-canonical form of the ML integrator. [29, 28] Note, 
however, that the system in Eqs. (42) is implicit in both X and Y. 

For the complementary generating function G(X,y) we follow the same procedure, 
defining ^(x,y) and ri(x,y) exactly as in Eqs. (28) and (29). We can now write Eq. (27) 
as 

d^(x,y)Ady = dXAdr](X,Y), (46) 
which, together with the definitions P = rj (X, Y ) and q = § (x,y), implies that 

G)' = -qdy-PdX (47) 
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is exact (dco' = 0). Therefore ©' = dG{X,y) = {dG(X,y)/dX)dX + (dG(X,y)/dy)dy and 

{(*,)~^, ,(i,n-^. («) 

The identity is given by Go(AT,y), where 

S(x,y) = - d -^, n( X ,y) = - d -^, (49) 
with -d 2 G /dxdy =p{x,y). We try G(X,y) = G (X,y)+hH(X,y), which leads to 



Z(x,y)=${X,y)-hu(X,y), rj{X,Y) = rj(X,y) +hv(X,y). 



(50) 



As in the case with Eqs. (42), this is the form we use as an integration scheme, and it is not 
necessary to obtain Go explicitly. Finally, a first order expansion in h again gives Eqs. (45), 
showing that Eqs. (50) are also a first order accurate integrator for Eqs. (20). 

The first order maps (x,y) — > (X, Y) defined by Eqs. (42) and Eqs. (50) involve integration 
of a one-degree of freedom Hamiltonian system with no explicit time dependence, and thus do 
not exhibit a parametric instability for small enough time steps. 

We have shown that both schemes in Eq. (42) and in Eq. (50) are first order accurate 
integrators for Eq. (20). What remains is to show that if they are applied sequentially one 
obtains second order accuracy. The first point is that if the maps based on F(x,Y) and G(X,y) 
are written as Tp (hi) and TgQi), respectively, one can show that 

T F l (h) = T G (-h) (51) 

and of course Tq (h) = Tp(—h). This follows from inspection of Eqs. (42) and (50). 
The composed map is T(h) — TqQi/2) o Tp(h/2) and we find 

T(h)- 1 = T F (h/2y l oT G (h/2)- 1 = T G (-h/2)oT F {-h/2) = T(-h). (52) 

Thus, T(h) is reversible and is therefore second order accurate. t We call T(h) = T G (h/2) o 
Tp(h/2) non-canonical symmetrized leapfrog (NSL). 

The methods outlined in this section preserve the non-canonical bracket, as does the actual 
time evolution of the non-canonical system. Integrators with this property are often called 
Poisson integrators[lS, 19, 20, 21]. 

It is interesting to note in this context that (fixed step size) CN can be written as 
TcN(h) — Tf(h/2) o Tf,(h/2), where Tf is forward (explicit) Euler and Tf, is backward (implicit) 
Euler. Further, the implicit midpoint method x' — x + (h/2) u (x) + (h/2)u(x! ) can be written 
as T mid (h) = T b (h/2) o 7/ (A/2). This shows that T mid = T f (h/2)~ l o T CN (h) o T f (h/2), so 
that the implicit midpoint method is related to CN by a non-canonical change of variables. 
This implies that the implicit midpoint method is equivalent to a symplectic method under a 
(non-symplectic) change of variables, so it is a Poisson integrator. Therefore it can be used to 
integrate a canonical Hamiltonian system with all the advantages of a symplectic integrator. 



3.4. Numerical results: non-canonical symmetrized leapfrog 

In this section we report numerical results where the non-canonical symmetrized leapfrog (NSL) 
integration scheme is applied to the cubic Hamiltonian [Eq. (18)], with p(x,y) = ft/Ao(l + 
a\x + a2y), as given inEq. (3.2). Integrating p(x,y) yields T](x,y) — hln(l+a\x + a2y)/Aoa2 

t Straightforward backward error analysis shows that, if a single-step map T(h) is accurate to first order and 
time-symmetric [i.e. satisfies T(h)~ i = T(—h)], then it must be accurate to at least second order. 
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and %(x,y) = h\n(l + a\x + a2y)/Aoai, and thus Eq. (42) with this Hamiltonian and time step 
variation becomes 

ln(l +a\x + a.2y') = ln(l -\-a\x-\-a2y) + Aoa%v{x), (53) 

ln(l -\-a\x +a2y') = ln(l -{-aix + aiy ) + AoaiM(;y), 

where u — dH/dy — y and v = —dH/dx = — (x + x 2 ). The first of these can be solved for 
3/ and then the second for x 1 , giving a non-canonical leapfrog integrator. The second non- 
canonical generating function method, given in Eq. (50), can similarly be written in terms of 
logarithms. As described in the last section, the composition of these two methods results in a 
symmetrized integrator, which can be written as: 

{l+a 1 x + a 2 y)e A ° a2vix y 2 - 



y = 



1 

x = 



■CL\X — 1 



/ai, 



(54) 



(1 +aix + a2y*) 



■a 2 y 



-1 



y 



1 +a\x +a2y )e' 



*\ A a 2 v(x')/2 



- a\x — 1 



fa\, 
/«2, 



Figure 7 shows numerical results from integrating the cubic oscillator using the NSL method 
of Eqs. (54) (shown in black) and fixed step size symmetrized LF (shown in gray). As for the 
extended phase space method, we find that the NSL method with this p{x,y) gives results with 
bounded errors in the energy. Again, this choice of p(x,y) was chosen because it suggests 
m = 1 resonance and leads to tractible equations; in Sec. 4.2 we show results with an optimized 
p(x,y) that minimizes the time-stepping error. The surface of section for this method looks 
qualitatively like Fig. 6, again with no resonant islands. 



o 

5f 
1 




Figure 7. The NSL method applied to the cubic oscillator, with time step variations as given 
in Eq. (3.2). Shown in black is the error in energy. For comparison, in gray is shown the 
error in energy for fixed step size symmetrized LF. Parameters: ho = 0.05, Ao = 1.0524, 
(90, po) = (0.3,0), a\ = 0.5, a 2 = 0.25. 



4. Using error estimates to find A(q,p) 

We next turn our attention to the choice of step size variation A = A(q,p), based on an error 
estimator. Basing the choice of A on an error estimate will give step size variations which are 
optimized to reduce errors, in a manner to be defined. 

4.1. Error estimates and minimization 

As an example of an error estimate we consider CN integration of dz/dt — u(z), where again 
z = (x,y). Namely, we define a map z' = 7},(z) given by 

z' = z + Au((z + z')/2), (55) 
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where z = z(f ) and z' = z(f + A), where A — At can vary as a function of z. If u comes 
from a Hamiltonian system, CN is symplectic, but that property is not important for these 
considerations. We can thus let A(f ) = A(q re f(t) , p re f(t)), where q re f and p re f are the exact 
orbits. fAs mentioned in Sec. 2.1, there will be terms proportional to dA/dt, but they will 
appear at the next order in A. 

The error estimate is found by expanding z(f + A) in A, and performing backward error 
analysis to obtain a local error estimate. This is a very different objective than that in the 
backwards error analysis in Sec. 2. In index notation we obtain 

dzi_ Adhi_ t£_d\i__ ( Adz A^(fz\ 

dt + 2 dt 2 + 6 *3 - u '{ z+ 2 dt + 4 dt 2 ) { ) 

A A 2 A 2 

= Ui{z) + -ZjUi.j + —ZjU L j + —ZjZ k u^ k . 



This leads to 



— = Uj + A' I —u k u jtk Uij - —UjU k u iJk I , (57) 
which can be written in vector notation as 

g=u + A 2 (^(u-Vu)-Vu-^u- ( Wu) ■ u) , (58) 
where u = u(z(f)). 

Writing the error per time step as e = A 2 w(t ), where w(t) — |(u ■ Vu) • Vu/12 — u • (Wu) ■ u 
the total error in an integration is 



E= edt= A z w(t)dt. (59) 
Jo Jo 

Now if we write the variable time step At as A = (dt / dx)h, where h = At is the constant step 
size in T, this error takes the form 

E = h 2 J dx(^j w{t). (60) 

We proceed to minimize the total error E by writing E — h 2 J dxJ£(t,dt/dx). Ignoring 
the constant factor h 2 , the Lagrangian Jz? is given by 

^{t,dt/dx)^{dt/dx) i w{t), (61) 

where z is the time, and we wish to find f (t) that minimizes E. This is done by forming the 
Euler-Lagrange equations, or equivalently by doing a Legendre transform to a Hamiltonian: 

dse „/ ' dt\ 2 . . 
P= dWdr) =3 \dr) (62) 

Jf(t,p)=pf^j-^f(t,dt/d'u). (63) 

We obtain Jff = 2(p/3) i / 2 w(t ) -1 / 2 . Notice that the Hamiltonian Jff is independent of the 
time T and is therefore a constant of motion. Also, the Hamiltonian J/f is proportional to the 

t Also, since we are considering A = A((), we do not need the terms proportional to H + po in the extended phase 
space method, Eq. (24). Regardless of the form of A, this term is of higher order than the error we are estimating in 
this calculation. 
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Lagrangian Jz? = (dt /dT) 3 w(t) = (p/3) 3//2 w(?) l l 2 , implying that Jz? is also a constant. Since 
dt/dx — 1 /p, we obtain an expression for the optimized time steps: 

p opt {z{t))=Cw{t) l l\ (64) 

This condition says that the Lagrangian, the integrand in Eq. (60), is constant with respect to T. 
(But the local error e in Eq. (59) is not constant with respect to t.) A condition of this form is 
called an equidistribution principle[3l]. 



4.2. Numerical results 



In order to test the optimized time step, we compare numerical results using constant time step 
[A = Ao], the optimized time step 

Kpt = C op ,w(q,py l/3 oc l/popt, 

and time steps designed to give equal arc length: 



(65) 



dH 



dH 
dq 



-1/2 



(66) 



With A arc , the magnitude of the phase space velocity v= (dq/dT,dp/dz) is constant. 
Numerical results for these three choices of A are shown in Fig. 8, using the CN integrator in 
extended phase space. For each of these methods, fixed steps in T of size At — h were used. 
The constant of proportionality C in A = C/p was then chosen so that, after /V periods of the 
orbit, t(N) equals t(N). This is done in order to distinguish the effect of varying the step size 
as a function of (q,p) from the effect of reducing the overall average step size. 

Two diagnostics of the numerical integrations related to the local error are shown in Fig. 8, 
for A = Ao (light gray), A arc (dark gray), and A opt (black). They are the local error A 2 w(f ) as 
a function of time and A 3 w(t), proportional to the integrand in Eq. (60), as a function of T, 
respectively. Variations in step size given by A arc reduce the error compared to the fixed step 
size integration, at almost all points along the orbit. Variations given by A opt further reduce 
the local error in some places, while allowing the local error to increase in others. As seen in 
Fig. 8b, this variation has the effect of spreading the local error along the orbit, so that A 3 w(f ) 
is constant as a function of T, as required by the equidistribution principle. 

As a further test to check if A op , really minimizes E [as defined in Eq. (60)], we also show 
numerical results where E is computed for 



An =Ca 



(1-P)+P*(9,P) 



-1/3 



[(l-j3)A + /3A o 



(67) 
(68) 

where Cp is chosen as above for each j3 , for a range of p 1 . When j3 = 0, the time steps are 
uniform, and when j3 = 1 the time steps are given by A upt . The results are shown in Fig. 9. As 
seen in the inset in Fig. 9, the error is indeed a minimum for j8 = 1, with E = 2.22 x 10~ 4 . For 
comparison, using A = A arc gives E — 2.48 x 10~ 4 , and A = Aq = const, gives E — 3.29 x 10~ 4 



5. Summary 

In this paper we have studied variable and adaptive time stepping for symplectic integrators, of 
importance for particle-in-cell codes, for accelerators, for tracing magnetic field lines, and for 
ray tracing. We have shown that problems observed in the literature for symplectic integrators 
with variable time steps fall into two categories. In the first, with time step A = A(f ), the 
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Figure 8. Comparison of the errors for the extended phase space CN method, for three time-step 
variations, A = const, (light gray), A arc (dark gray), and A op , (black), (a) The local error as a 
function of time t [integrand of Eq. (59)], and (b) the local error as a function of X [h 3 times 
the integTand of Eq. (60)]. Parameters: h = 0.1, (qo,Po) = (0.4,0), A = 0.1, C aTC = 0.41, 
C pt — 0.28 
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Figure 9. The global error [Eq. (60)], over the time T = 20, for the time step variation A^ 
[Eq. (67)]. As seen in the inset, the error is minimized with E = 2.22 x 10~ 4 for p = 1, for 
which Ap = A opl . For comparison, using A = A olr gives E = 2.48 x 10~ 4 , and A = Ao gives 
£ = 3.29x 10" 4 . 



integrators are still symplectic but results can exhibit parametric instabilities associated with 
resonances between the time step variation and the orbital motion. We have characterized these 
instabilities by means of backward error analysis and numerical integrations. In the second 
category, the time steps depend on position in phase space, A = A(q,p), and integrators which 
are symplectic for A = const, are no longer symplectic. 

We have discussed symplectic integrators for A = A(q,p) of two basic varieties. In 
both, we start by modifying the time variable t X such that p(q,p)dt — dx and p(q,p) 
A(q,p) and we do uniform time stepping in x. In the first method, one extends the phase 
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space (q u ...,q n ,p u ...,p„) — ► (q ,qi,...,q n ,p ,pi,...,p n ), with q = t and p its canonical 
conjugate. In this extended phase space with t taking the place of time and At = const., any 
numerical integrator can be used and we investigate modified leapfrog[28, 29], symmetrized 
to be of second order, and Crank-Nicolson. In the second variety of integrator, we start 
by noticing that the equations to be solved, written as functions of T, are in Hamiltonian 
form with non-canonical variables for one degree of freedom problems. We then construct 
a Poisson integrator (the analog of a symplectic integrator but for non-canonical variables) 
in terms of a generating function of mixed variables, which is a generalization of the mixed 
variable generating functions for canonical variables[32]. We also show how to symmetrize this 
integrator to obtain second order accuracy, obtaining the non-canonical symmetrized leapfrog 
integrator. 

We have investigated these integrators for a model Hamiltonian system (the cubic 
oscillator) and found no parametric instabilities or growth or decay due to lack of phase 
space conservation. 

We have also investigated adaptive time-stepping, focusing on minimizing the norm of the 
error obtained by backward error analysis. We have applied this optimal time-stepping to the 
model Hamiltonian system and found that indeed the error is minimized while the advantages 
of symplectic integration are preserved. 
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Appendix A. Variable time steps and perturbation theory 

Two main time step variations are considered in the literature: A = A(f ) and A = A(q,p). The 
phase space dependent variation A(q,p) is natural because the errors in an integration method 
of an autonomous system depend only on phase space position. Indeed, we found that error 
analysis (as in Sec. 4) for the integration method gives errors which depend on location in 
phase space (q,p). 

The second kind of variation A(f ) comes from a type of perturbation analysis. The logic 
is as follows: orbits (q(t),p(t)) of a Hamiltonian system are often periodic, and because of 
this, time step variations A(q,p) on the orbit should also be periodic. That is, we can substitute 
q(t) and p(t) to arrive at A(f) = A(q(t ),p(t)), giving periodic time dependence. 

This argument ignores the stability issue: if an orbit is perturbed, the period of A(q(t), p(t)) 
will vary, but this will not occur for a prescribed A(t ). 

It has been observed[5, 7, 8, 9] that such a A(f ) is not a stable method for varying the 
time step of a symplectic integrator. Some analysis (see Ref. [16] and Sec. 2) shows that this 
does not work because it introduces parametric resonances into the system. As we have shown 
in Sec. 2, the integrator is still symplectic in this case, but the problem is that the modified 
equations (in the case of a harmonic oscillator) take the form of the Mathieu equation, which is 
still Hamiltonian but can be unstable. 
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We here review how such parametric instabilities and secularities arise in perturbation 
theory, to help us understand the parametric instabilities of Sec. 2.1. 

Secularities show up in perturbation theories in the following manner. Take the system 

x+x + ex?=Q. (A.l) 



A straightforward perturbation method involves iterating 

h+\ +x, 

For xq = bcost, we obtain 



h+\ +Xk+\ = -ex\. (A.2) 



ii +x\ — — eb 3 l~cos3f + ^cosfj . (A.3) 

The last term on the right is resonant and leads to a secularity, 

3sb 3 

x\ = bcost — fsinH (A. 4) 

8 

However, it is clear from Eq. (A.l) that the exact solution is bounded: the energy is 
H = x 2 /2 + x 2 /2 + ex 4 1 4. The t sin t term represents the frequency shift that occurs in Eq. ( A. 1 ) 
because of the term proportional to £, and is accurate for short time, but clearly wrong for long 
time. 

For an example of a perturbation approach with a parametric instability, we consider 
instead 

x+x + ex 2 = 0, (A.5) 
with the perturbation method 

Xk+i +Xk+\ = sx k x k+ i . (A.6) 

(We choose a different model in order to obtain again results at first order in the perturbation 
theory.) From xq = bcost we obtain 

ii + (l+efecosf)jci =0. (A.7) 

Using T = t/2, this is the Mathieu equation d 2 x\/di 2 + (a + 2qcos2t)x\ = 0, with a = 4 
and q = 2eb, for which a parametric instability occurs. As in the secular example above, the 
exponential growth represents the correct frequency shift for small time, but is wrong for longer 
time. Indeed, the energy is H = x 2 /2 +x 2 /2 + ex 3 /3, and for initial conditions near x =x = 
the solution cannot grow indefinitely. As we have seen in Sec. 2.1, parametric instabilities for 
Af = A(f ) have a very similar character. 

Appendix B. Non-canonical brackets and forms 

In this appendix we first describe the conditions that change of variables must satisfy so that 
the form of a set of equations defined with a non-canonical bracket is preserved. We then show 
that the flow of a Hamiltonian system written in non-canonical variables, such as Eq. (20) (one 
degree of freedom), generates a change of variables satisfying these conditions. 
First define the non-canonical bracket 

r , x , v, 1 df(z) dg(z) 

In one degree of freedom, [•, ] z is indeed a non-canonical bracket, which can be verified by 
proving the Jacobi identity. In more that one degree of freedom, this does not define a legitimate 
non-canonical bracket, since it does not satisfy the Jacobi identity. 



1 _ 3 
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We now ask what change of variables Z = Z(z) preserves the form, i.e., what conditions 
must this change of variables satisfy so that 

[/(z),g(z)] z =[/(Z),g(Z)] z , (B.2) 

where 

p(ZJ dZ m dZ n 
Using the definitions of these brackets and the chain rule, Eq. (B.2) becomes 

1 Ct dZ m df(x) dZ n dg(z) = 1 ^ df(Z)dg(Z) (R4) 
p(z) '■' dzi dZ m dzj dZ„ p(Z) '"" dZ m dZ n 
In order for this to be true, we require that 
1 dZ m dZ„ _ 1 



MEM 1 = ^j^E, where A/y = E tJ = £, 7 . (B.6) 



p(z) u d Zi dzj p(Z) 
or, written in matrix form, that 

t P ( z ) dZj 

= WrT\ E > where M '> = IT 
p(Z) ■' dzj 

If this condition is satisfied, then the transformation z — > Z preserves this specific non-canonical 
symplectic structure. 

The next issue is to show that the time evolution operator for Eq. (20) generates a map 
with this property. Since Eq. (20) can be written in terms of the non-canonical bracket as 

g = (B.7) 
then, for infinitesimal h, if we define zi = Z;(0) and Z, = Zi[h) we have 

Zi=Zi + h[zi,H] z . (B.8) 
Taking the bracket of Z, and Zj therefore gives 

[Zi,Zj] z = [zi,Zj} z +h[[zi,H} Z) zj] z +h[zi,[zj,H] z } z (B.9) 

= [zi,Zj] z +h[[zi,H} z ,zj] z + h[[H,zj] z ,Zi\ z (B.10) 

= [zi,Zj] z +h[[zi,zj} z ,H} z . (B.ll) 

The last equality follows from the Jacobi identity. This equals 

1 



[Z h Zj] z = Bij ( -^+h[l/p,H] z ) (B.12) 



-^ + (Z k -z k )^- (-*--)) (B.13) 
p(z) dzk\p{z)JJ 

^ W) = [Z - Z]]z - (B - 14) 

Thus, the time evolution operator for an infinitesimal time step h preserves the non-canonical 
symplectic structure. This fact is the basis for finding a non-canonical generating function in 
Sec. 3.3 that gives a first order accurate integrator for the non-canonical Hamiltonian system. 
Finally, we express these ideas in terms of differential forms. Eq. (B.5) can be written 

, . dZi- dZi , . 
P (Z)^e kl ^=p(z) £ij . (B.15) 
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From this, we obtain the expression 

p(Z)(dZ k /dzi)e k ,{dZ,/dzj)dzidzj = p(i.)eijdzidzj, (B.16) 

or 

p (Z)e M dZ k dZ l = p {z) £ij dzidzj. (B. 17) 

Since we are considering only 2D phase space, this can be expressed as 

p{X,Y)dXAdY =p(x,y)dxAdy. (B.18) 
Eq. (B.18) is the two-form used in Sec. 3.3. 
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